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Abstract. Many sociological networks, as well as biological and technological ones, can be represented 
in terms of complex networks with a heterogeneous connectivity pattern. Dynamical processes taking 
place on top of them can be very much influenced by this topological fact. In this paper we consider a 
paradigmatic model of non-equilibrium dynamics, namely the forest fire model, whose relevance lies in 
its capacity to represent several epidemic processes in a general parametrization. We study the behavior 
of this model in complex networks by developing the corresponding heterogeneous mean-field theory and 
solving it in its steady state. We provide exact and approximate expressions for homogeneous networks 
and several instances of heterogeneous networks. A comparison of our analytical results with extensive 
numerical simulations allows to draw the region of the parameter space in which heterogeneous mean-field 
theory provides an accurate description of the dynamics, and enlights the limits of validity of the mean-field 
theory in situations where dynamical correlations become important. 
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1 Introduction 

The heterogeneous topology of a networked substrate has 
been proven to have a large impact on dynamical pro- 
cesses taking place on top of it |l|2j . These topological 
effects are especially remarkable in the case of scale-free 
(SF) complex networks |3|4j , characterized by a degree dis- 
tribution P(k), defined as the probability that an element 
in the network (vertex) is connected to k other elements, 
that exhibits a power-law behavior, P(k) ~ fc -2-7 , with 
< 7 < 1 Q. The diverging second moment (fc 2 ) of the 
degree distribution has thus been found to be at the core 
of the peculiar behavior observed in a wide array of non- 
equilibrium dynamical processes, ranging from percolation 
|5l6j . absorbing-state phase transitions [7], self-organized 
criticality |8|9j . synchronization phenomena 10 . opinion 
dynamics [TT], etc. 

The interplay between topology and dynamics has been 
particularly studied in the case of epidemic processes [12] , 
where the relevant substrate is the network of contacts 
through which the disease spreads [13] . Starting from the 
first observations of an epidemic threshold scaling as the 
inverse of the second moment of the degree distribution, 



1 To ease the notation in our mathematical treatment, we 
will use this definition of the degree exponent in a power-law 
degree distribution. 



and thus vanishing in the thermodynamic limit of an in- 
finite network size |14I15I16I17| . a wealth of interesting 
and relevant results have arisen, dealing, to mention just 
a few, with immunization strategies |18ll9j , effects of bi- 
partite (heterosexual) populations |20| or epidemic fore- 
casting [2"T] . 

The understanding of the features of epidemic spread- 
ing is mainly based on the analysis of compartmental mod- 
els |22j . in which the population is divided into different 
classes, according to the stage of the disease. Individuals 
(the vertices in the network) are in this way classified as 
susceptible (healthy and capable to contract the disease), 
infected (sick, and capable to transmit the disease), re- 
covered (immunized or dead), etc. With these definitions, 
different epidemic models can be formulated, according to 
the succession of states that the evolution of the disease 
imposes on each individual, such as susceptiblc-infccted- 
susceptiblc (SIS), susceptible-infected- recovered (SIR), su- 
sceptible-infected-removed-susceptible (SIRS), etc. The the- 
oretical analysis of the behavior of these compartmental 
models in complex networks starts from the application 
of the heterogeneous mean-field (HMF) theory [112] . This 
formalism is based on the assumptions that all vertices 
with the same number of connections (i.e. within the same 
degree class) share the same dynamical properties, and 
that fluctuations are not important, and therefore all rel- 
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evant variables can be described in terms of deterministic 
rate equations. The first assumption becomes natural once 
we admit that the degree is the only parameter describing 
the state of a vertex. On the other hand, the second as- 
sumption finds support in the small-world property shown 
by most complex networks |23| , implying that dynamical 
fluctuations take place so close together that they can be 
washed away in very few time stepfl HMF has proved 
to be extremely useful in providing an accurate descrip- 
tion of epidemic models on complex networks, and has in 
fact become the de facto standard tool to analyze general 
non-equilibrium processes on such substrates p] . 

In this paper we will pay attention on a non-equilibrium 
dynamical model with relevance both in epidemic model- 
ing and other ambits of non-equilibrium statistical physics, 
namely the forest fire model (FFM). First introduced in 
1992 by Bak et al. [25], and further developed by Drosscl 
and Schwabl [26], the FFM was elaborated to show self- 
organized criticality and avalanche behavior in a specific 
limit of its defining parameters. Even though its general 
status as a self-organized critical model is under debate 
[27128] , it has found successful applications as a general 
disease propagation model [12129] . in which susceptible 
individuals can get the disease either by transmission from 
an infected neighbor or spontaneously (because of a mosquito 
bite for instance) , while recovered individuals can become 
again susceptible. It is thus akin to a SIRS model [T2] 
with an external source of infected individuals. More in- 
terestingly, it encompasses several other compartmcntal 
models, which can be recovered in a convenient way as 
certain limits of the parameters that define the FFM. 

Previous works on the FFM in complex networks have 
reported, among other results, the presence of self-sustained 
oscillations in small- world networks [30] and analyzed the 
distribution of excitations depending on topology [31] . From 
the perspective of the SIRS model, on the other hand, its 
epidemiological implications have been discussed for cer- 
tain ranges of its parameter space [32133] . In the present 
paper we provide an extensive theoretical analysis of the 
FFM, using the HMF formalism and focusing on the steady- 
state dynamics of the model in the whole range of its pa- 
rameter space. Our analysis allows to emphasize its in- 
terpretation, in the different regimes, in terms of known 
epidemic models, providing analytic expressions for the 
steady-state density of infected individuals in certain lim- 
its of the relevant parameters. A comparison of the theo- 
retical results with extensive numerical simulations, allows 
finally to unveil the limitations of the HMF approach in 
this and probably other epidemic models, hinting towards 
the break down of HMF theory when dynamical correla- 
tions become relevant [24] , 

We organized our paper as follows: In Sec. [2] we define 
the FFM, discussing its relation to self-organized critical- 
ity and disease propagation. In Sec. [3] we develop the HMF 
theory of the FFM in general complex networks. Sec. [4] 
deals with the steady-state solution of the HFM equa- 
tions obtained before. We show in particular how an ap- 




h (n.n.) 



Fig. 1. Representation of the three states present in the for- 
est fire model, and of its dynamical rules, (n.n. stands for a 
"nearest neighbor" interaction). 



propriate rescaling of the equations allows to simplify the 
description and to conveniently reduce the number of de- 
grees of freedom. A general analysis is presented for both 
homogeneous networks and heterogeneous networks with 
no spontaneous infection. An explicit analysis of uncorre- 
cted SF networks is presented in Sec. [5] The numerical 
simulations shown in Sec. [6] allow us to check the validity 
of our theoretical results, as well as to draw the limits of 
validity of general HMF approaches. Finally, we present 
our conclusions in Sec. [7] 



2 As a matter of fact, fluctuations can be shown to be irrel- 
evant in some particular cases [24] , 



2 Forest fire model on complex networks 

We consider the FFM on general complex networks which, 
from a statistical point of view, are described at a coarse- 
grained level by the degree distribution P{k) and the degree- 
degree correlations, given by the conditional probability 
P(k'\k) that a vertex of degree k is connected to a vertex 
of degree k 1 [34] . 

In the FFM, each vertex in the network is in one of 
three excluding states: E (empty), T (tree), F (burning 
tree) . The evolution of the model is defined in a continuous 
time formulation in terms of the possible events that can 
happen in a small time interval At (see Fig. [IJ 

1. E — > T : A tree can grow on an empty vertex with 
probability pAt. 

2. T — > F : For a tree, each of its burning neighbors (if 
any) can light it with probability hAt. 

3. T — > F : Additionally, there is also a probability gAt 
for a tree to catch fire spontaneously (e.g. mediating 
a lightning). These two burning events are considered 
probabilistically independent. 

4. F — > E : A burning tree leaves an empty vertex with 
probability I At. 

The FFM was proposed to exhibit self-organized crit- 
icality in the double limit p — ► and g/p — > 0, in which 
clusters of trees are allowed to grow before burning down, 
leading to a distribution of fire avalanches with an approx- 
imately power-law form [35] , For finite p and g values, an 
activated dynamics with no evidence of avalanches is in- 
stead observed |30|31] , 
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From the point of view of disease modeling, trees, fires, 
and empty sites represent, respectively, susceptible, in- 
fected, and recovered individuals. Within this interpreta- 
tion of the model, susceptible individuals can acquire the 
disease by contact with one or more infected individuals 
at rate h; infected individuals recover in a time scale of the 
order ~ 1 /I; recovered individuals become again suscepti- 
ble in a time scale ~ 1 /p; and healthy individuals become 
spontaneously infected in a time scale ~ !/<?. With a set 
of 4 parameters, a unified description of the most usual 
epidemic models is achieved. Thus, in the limit p = g = 0, 
we recover the SIR model; the limit g = and p — > oo 
leads to the SIS model; while the limit g = corresponds 
to the SIRS model. 

Note finally that our definition of the FFM considers 
time as a continuous variable, in opposition to previous 
approaches. This formulation is preferred in order to lead 
more naturally to a continuous analytical description in 
terms of differential equations, and will be taken into ac- 
count when performing numerical simulations in Sec. [6] 



3 Heterogeneous Mean- Field theory for the 
FFM in complex networks 

Within the HMF approach, a dynamical system is as- 
sumed to be fully determined in terms of the relative prob- 
abilities that a vertex of given degree k is in any one of the 
states allowed by the dynamics |14ll6j . In the case of the 
FFM, this description invokes the partial densities p%(t), 
defined as the conditional probability that a vertex of de- 
gree k is, at time t, in the state a, with a e {E,T,F}. 
Since each vertex must be in one of these states, the par- 
tial densities satisfy the normalization condition 



Pu(t)+p T k {t) + p F k (t) = l. 



(1) 



Therefore, only two independent partial densities, say pf (t) 
and Pk(t), must be considered in the analysis. On the 
other hand, the density of vertices in each state at time t 
is given by 

P a {t) = Y. p ^)p a ^ t )- ( 2 ) 

k 

At the core of the HMF theory lie the rate equations 
fulfilled by the partial densities. By considering the differ- 
ent microscopic steps allowed in the model we can readily 
write the change of the quantities p k (t) in an infinitesimal 
time step At, that are given by 

Pfc (t + *t) = Pk(t) + Pk (t)Prob k (F - E) 

- pi (i)Prob fe (£ -+ T), (3) 
p T k (t + At) = P l(t) + p*(t)PToh k (E -> T) 

- pl(t)Prob k (T ^ F), (4) 
pi (t + At) = pi (t) - pi (t)Prob k (E -> F) 

+ pi(t)Pvoh k (T^F), (5) 

where Probfe(a — + (3) is the probability that a vertex of 
degree k experiences the transition from the state a to 



the state p in a time interval At. From the definition 
of the FFM in Sec. [21 we can immediately write down 
Prob^E 1 — » T) = pAt and Prob fc (F -> E) = I At. In 
order to construct the term Probfc(T — > F), we must 
consider that rules 2 and 3 defining the model, which 
represent a tree catching fire, are statistically indepen- 
dent. Therefore if we define H as the event "A tree is 
lighted by its neighbors" and G as the event "A tree 
lights up spontaneously" , we have that Probfc (T — * F) = 
Prob fc (GU#) = Prob fe (#)[l-Prob fc (G)]+Prob fe ((7). Rule 
3 gives Probfc(G) = gAt. Now, since in this description 
vertices with the same degree are statistically equivalent, 
the state of a given vertex is independent on the state of 
its neighbors, and it only depends on its degree. This al- 
lows us to write the probability for one neighbor of a tree 
with degree k to be burning as 



9 k = J2P(k'\k)pi, 



(6) 



given in terms of the average of the conditional probability 
P{k'\k) that the vertex k is connected to a vertex of de- 
gree k' , times the probability that this last vertex is burn- 
ing, pi . Notice that here we are assuming that the edge 
through which k 1 became burning is immediately available 
to transmit again the fire. This assumption, in opposition 
to the behavior of the SIR model [36] , will thus be valid 
only for p > 0. 

From rule 2, the probability that a particular nearest 
neighbor fire ignites a tree in a vertex of degree k is given 
by hAt9 k . Therefore, the probability that a tree of degree 
k is ignited by any of its nearest neighbors is PTob k (H) = 
1 - [1 - hAt6 k ] k . Substituting this expressions in Eqs. ([3])- 



([5]) , and taking the limit At 
equations for the FFM, 



0, we obtain the final HMF 



Pi(t)=Pi(t)~P P i(t) 

pi(t)= P pi(t)-pi(t) [hkj: k ,p(k'\k)pi(t)+ g } , 

Pi(t) = PK*) [hkEk' P(k'\k)pi(t) + g] pi(t) 

(7) 

where we have set I = 1, which amounts to a trivial rescal- 
ing of time. 

Eqs. ([7]) represent a complete description of the FFM 
at the HMF level. Even though we have derived them in a 
phenomenological way |I4|37j . they can also be obtained 
from a microscopic point of view, considering explicitly the 
state of each vertex evolving as a Poisson random process 
|24|38j . The mean-field result is then recovered by averag- 
ing over the random processes and over the vertices with 
same degree. 

One final warning comment is in order here, concerning 
the fact that, in writing Eqs. ([7]), we have neglected alto- 
gether dynamical correlations between adjacent vertices, 
assuming explicitly that the state of a vertex is indepen- 
dent of the state of its nearest neighbors. As we will see, 
this assumption is not correct, especially for low fire (in- 
fection) densities, when the positions of different fires are 
in fact strongly correlated, leading thus to a breakdown 
of the HMF theory predictions (see Sec.[6|). 



4 



J.-D. Bancal and R. Pastor-Satorras: Steady-State Dynamics of the FFM on Complex Networks 



©gb© 

h (n.n.) 

Fig. 2. Representation of the two states present in the SIS+g 
model, and of its dynamical rules. Its steady-state is directly 
related to that of the FFM. (n.n. stands for a "nearest neigh- 
bor" interaction). 

4 Steady-state solution in general networks 

Let us consider the long time properties of the FFM. It 
typically corresponds to the steady-state calculated by set- 
ting Pk = Pk = Pk = Vfc in tnc HMF E( l s - ©■ which 
yields the algebraic equations 



F E 
Pk -PPk 

PP E k-P T k [hkEk'P(k'\k)p£ 



(8) 



Pi [hkEk'P(k'\k)pF,+g] 



Pi 



for the (now time- independent) variables p%. We look for 
nontrivial steady-states, so we will be concerned in search- 
ing solutions with p% ^ 0. 

The analysis of Eqs. |(HJ) can be simplified by noticing 
that the empty state plays the role of a rest state (c.f. Ap- 
pendix A), which the system enters and leaves with con- 
stant rates. It can thus be factorized by writing its popu- 
lation density, from the first equation in ([8]), as p F = p F jp 
and substituting it in the other two equations. Therefore, 
introducing the rescaling factor 



77=1 



1 

P 



and defining the new variable and parameter 

9 = TO. 



Pk = VPk, 



we can consider the simplified set of equations 



Pl-Pl [hkEk'P(k'\k)pF 
Pk+Pl 



(9) 



(10) 



fill 



as characterizing the steady-state of the FFM in a general 
complex network with a correlation pattern given by the 
conditional probability P(k'\k). 

From Eq. (jTTJ) , we can see that the steady-state of 
the FFM can be in general mapped to the steady-state 
of an SIS model [2] with a random source of infected 
individuals, arising from isolated susceptibles with rate 
g. We can thus refer to it as a SIS+g model, see Fig. O 
In particular, setting 5 = 0, the FFM becomes the SIRS 
model, which is therefore exactly mappable to the SIS 
model [32]. 

In order to solve the set of Eqs. (fll|) , one can proceed to 
substitute its second equation into its first one, to obtain 

Pi 



F as a function of 9k = i]@k, namely 



Pi 



hkOi 



9 



1 + {hkO k + g) 



(12) 



The equation is closed by expressing 9k self-consistently 
as 



=yp(k'\k) ru f Sv . +3 



(13) 



Solving this equation for 9k directly gives the steady-state 
fire density, and thus the stationary solution of the pro- 
cess. 

Prior to solving these equations, however, notice that 
Eqs. (|11|) imply that the solution for p F takes the func- 
tional form p F (g, h), and is now independent of p (namely 
of if). Therefore, we can write down a scaling solution for 
the steady-state fire density in the FFM in any network 
as 

p F ( V ,g,h) = U F (r,g,h). (14) 

This scaling solution implies that the factor p (growth of 
new trees) affects the model only by a rescaling of the 
fire density and of the rate of spontaneous lightning. This 
justifies the fact that no more than 2 out of the 4 ini- 
tial parameters are relevant to the study of the stationary 
state of the FFM within the HMF theory approximation. 
For the particular case g = 0, Eq. ([14")) leads to 



p F (r,,g = 0,h) = -F(h), 



(15) 



that is, the fire density is a function of h, divided by r\. 
In the limit p — > oo (77 — > 1), we recover the standard 
SIS model. On the other hand, any finite p will lead to a 
smaller value of the fire density. The fact that an upper 
bound on p F can be deduced from 77 is in general valid for 
any g: since p F < 1, given that it is a probability, we have 
that 



1 

< - = 

T) 



P 



p+1 



(16) 



This upper bound on p F has important consequences from 
a numerical point of view. In fact, even in the regions of 
the parameter space (h, g) where a nonzero value of p F is 
expected, a very small value of p will lead to a correspond- 
ingly small fire density, which can be difficult to measure 
unless in the limit of very large network sizes. 

The factorization of the empty state, the functional 
form of the fire density and its upper bound in terms of 77 
are in fact quite general features, that can be found in any 
dynamical system having rest states, see Appendix A. 



4.1 Homogeneous networks 

Let us consider first the simplest situation of the FFM 
taking place in a homogeneous network, in which the de- 
gree distribution is peaked at an average degree (k) and 
decays exponentially fast for k <C (k) and k ^> (k). We 
can approximate all vertices as having the same degree 

p F effectively inde- 



k = (k). In this case, we have p F 



pendent of k, and also 9k 
form 



p b . Eq. ([T2| thus takes the 



F\1 



h(k)(p*) 



1 



h(k))p* 



0. 



(17) 
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whose only positive solution is 



- F h(k) 
P = 



l-g + y/-4h(k) + (l+g + h(kjy 
2h(k) 



(18) 



as already shown in Ref. [39] . 

For g > the fire density is strictly positive. This 
results from the fact that the trees in the network can 
always ignite themselves with some nonzero probability, 
therefore always reviving the fire density. On the other 
hand, for g = 0, the fire density takes the form 



|l-fe<fc)|-(l-fe<fc» 
2r)h(k) 



(19) 



which is equal to for h < and positive otherwise. 

That is, the FFM experiences an absorbing-state phase 
transition [40] at a critical value h c = (k)^ 1 . For h > h c , 
the systems is in an active phase, in which the fire activity 
never stops, taking the asymptotic form 



h — h c 



On the other hand, for h < h Ci the system reaches an 
absorbing state, in which fire always ends up disappearing 
by lack of transmissibility. 



From Eq. (|12[) . is now an algebraic function of k. For 
the interesting case of SF networks, with a degree distri- 
bution P{k) = (7 + l)r7i 7+1 fc~ 2 ~ 7 in the continuous degree 
approximation, where m is the minimum degree present 
in the network, Eq. (|13[) reads, replacing summations by 
integrals, 



7m ' 



h6k-~< + gk 



-1-7 



F( 1.7J7+1 
75 

(7 + l)mh0 J 



1 + 7]{hk6 + g) 
mh6 ) 



F 1 )7 +1;7- 



1 + 5 
mhO 



(23) 
(24) 



where F (a, 6; c; z) is the Gauss hypergeometric function 
[12] • Using the power series development of the hypergeo- 
metric function [42] or the asymptotic expression 



(20) F(i, 7 ;i + 7;-^ 1 )=7E(- 1 )" — 



77T 



n=l 



7 sin(77r) 



(25) 

valid for arg(z) < 7r and 7 ^ N, one finds the self-consistent 
equation for 9 to be: 



4.2 Phase transition for g = on complex networks 

For networks with a general degree distribution P(k) and 
general correlation pattern P(k'\k), the explicit solution 
of Eqs. (fT2"[) and (flU)) becomes a quite difficult task. It is 
possible, however, to obtain information for a general net- 
work in the particular case g = 0. Setting g to changes 
the forest fire model to a SIRS model, whose stationary 
state can be related to the one of the SIS model by removal 
of the rest state, see Eq. (jTTj) . Therefore, in this particular 
limit, the FFM exhibits an absorbing-state phase transi- 
tion between an active (burning, infected) phase and an 
absorbing (fire-free, healthy) phase, located at the critical 
point [41] 

hc = ~l-, (21) 

where A m is the largest eigenvalue of the connectivity ma- 
trix Ckk' = kP(k'\k). Interestingly this threshold is inde- 
pendent of the rate p of creation of new trees, which only 
affects the overall density of fires, as expressed in Eq. (fTS"!) . 



9 



1 ( 1 

—F 1,7;7+1;- 



mh6 



(26) 



for all > 0. This is the final equation we need to solve 

mhO ^ 

in order to find the steady state fire density. Note that the 
condition of validity 7 ^ N is no more a restriction here, 
by analytical continuation of the hypergeometric function. 

Before proceeding, let us express more explicitly the 
dependence of the fire density (what we are ultimately 
interested in) and the probability for a neighbor to burn, 
9. We can directly calculate it using Eq. (fl"2"|) and a similar 
reasoning as before, to obtain: 



£P(fc)/f = ( 7 +lW +1 



hk9 + g 



1 



1 + 5 



F l,7+l; 7 + 2 



1 + r)(hk6 + g) 
1+5 



mhO 



(27) 



Now, using one of Gauss's relations for contiguous hyper- 
geometric functions [42] . namely 



5 Explicit solution for uncorrelated scale-free 
networks 

In order to obtain explicit analytical results for the HMF 
equations of the FFM, we restrict ourselves to the case of 
uncorrelated networks. In this case, the conditional prob- 
ability P(k'\k) takes the form P(k'\k) = k'P(k')/(k) 0]. 
9k thus becomes independent of k: 



= 9 = 



1 



Y,kP(k)p 



{a—c)zF(a, b; c+1; z)+cF(a, 6—1; c; z) = c{\—z)F(a, b; c; z) 

(28) 

with a = 1, 6 = 7+1 and c = 7 + 1, the hypergeomet- 
ric function can be re-expressed directly in terms of 9 as 
obtained in Eq. ([26]) to get 



fJ 



7 + 1 mh9 



1 



1 



1 



(1-0) 



(29) 



(22) ^ n a g enera l SF network without correlations, the fire den- 
sity is thus a quadratic function of the probability 9. 
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5.1 Exact solution for 7 = 1 

In the case 7 = 1 it is possible to solve exactly the self- 
consistent equation Eq. ([26]) . which takes the form [42] 



(30) 



1 + 3 (I+3) 2 



mh6 ( l + o 

In ( H 1 

mhO 



Introducing the new variable y = 3 + ^fjjr the equation 
becomes 

v _ g 1 (!+g) 2 



(31) 

Recognizing the solution of ye y = x as the W Lambert 
function, y = W(x) [43], the result follows: 



1 + 5 
mft 



-l + g- 1 W[ge 



(32) 



The value of the fire density is then obtained by plug- 
ging this expression into Eq. ([29]) . Expanding the fire den- 
sity at first order for small 3 yields then 



2hm I e * 



1 



hm (ehm. — lj 

8 - Whm 
h 2 m 2 [e*^ — 1 



hm 



8 - 2hm 
hm ( ew — 1 



(33) 

where the first term corresponds, obviously, to the SIS 
result |14j . recovered in the limit g — > 0. The expansion in 
terms of ft, is 



.9 



2gm 



1 + 5 (l + .g) 3 ' 



(34) 



yielding a nonzero fire density (an infected steady-state) 
for any value of ft. if 3 > 0. 



5.2 Asymptotic solution for 7 7^ 1 

Let us now turn our attention to the behavior of the FFM 
in the general case 77^ 1. To do so, we will study the limit 
of low fire density, namely 9 -C 1 and g -C 1. In this limit, 
Eq. ([29]) leads to 



7+1 

7 



mh9 1 



(35) 



so we need to develop the self-consistent equation for 9, 
Eq. (|26p . up to the first most relevant terms, using the 
asymptotic expansion Eq. ([23]) . For g = we recover the 
known result p F ~ (ft — ft c )^ [37], where 



ft r; = 



for < 7 < 1 
for 7 > 1 



(36) 



and 



Y+^ for < 7 < 1 
/? = for 1< 7 < 2 

1 for 7 > 2 



(37) 



For g > 0, let us consider separately each possible value 
of 7. 

(1) < 7 < 1: 

In this case, the leading approximation for 9 is 



77T 



mh9 



1+.9 (1 + 3) sin(77r) \1 + g 



(38) 



which is valid for <C 1. If 3 7^ 0, the solution of 9 
depends in a nontrivial way on both g and ft. as 



77T 



1 + g ' (I+3) sin(77r) 



3 



(39) 



where the function k(o, b) is defined as the solution of the 
implicit equation n(a, b) = a + b ■ «(a, b) 1 . In Appendix B 
we give an explicit expression for k(o, b) in terms of a and 
b. For small 3 and ft, two different regimes can be isolated 
from the development of the k function ([7T]) : 



— ft <C 3 T : The most significant terms are now 



9 



77T 



mhg 



(40) 



1 + 3 (1 + 3) sin( 7 7r) \(1 + g) 2 , 
so that, to leading order in ft, the fire density goes like: 
-f 9 , 7 + 1 m 9 



1 



7 



-ft. 



(41) 



The introduction of an infinitesimal g > destroys 
again the absorbing-state phase transition, since both 
terms in Eq. (|41| arc positive. 
3 <C ft 1 "^ : In this case, the leading terms are: 

00 + (1 - 7 )- 1 3, (42) 
whcrc ^0 = ( 2 ^gf ) ^- Thus, Eq. pj yields 



P 



7+1 

7 



mfte»°(l-6» ) 



1 



(1 + 7)mft 
7(1 - 7) 



(1 - 29 Q ) 



(43) 

One should notice that in order to derive this last ex- 
pression we used Eq. ([3"8]) . which assumes j^i <C 1, so 
not only do we need g to be small in order for this ap- 
proximation to hold, but ft cannot be too large either. 

(2) 7 > 1: 

The self-consistent equation Eq. ([2"o| can be approxi- 
mated in this regime by: 



9 



7 



ih9 



1 + 3 7-1(1+ gf 



(44) 
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So wc obtain: 



( 7 -l).9(l + .g) 



(7 — 1)(1 + g) 2 — jmh 
and to first order, the fire density is given by 
f ~g(l-g + 2mh). 



(45) 



(46) 



To summarize, in all cases considered above the pres- 
ence of a fire lightning probability g > eradicates any 
phase transition in the model, just like on homogeneous 
networks, Sec. 14.11 rendering a nonzero steady-state that 
grows, at lowest order, linearly with g, with a numerical 
prefactor which is a complex function of h, depending on 
the particular value of 7 considered. 
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6 Numerical Simulations 

In order to check the analytic predictions developed in the 
previous sections, we performed extensive numerical sim- 
ulations of the FFM on top of different network models, 
both homogeneous and heterogeneous. Simulations were 
implemented using a sequential update algorithm |40j . 
Given a substrate network, each vertex is first randomly 
initialized in one of the three possible states, empty, tree, 
or fire. In the dynamics, every time step a vertex is cho- 
sen at random and its state is updated applying the rules 
defined in Sec. [2] with given reaction rates £, p, g and h. 
In principle, to reproduce the exact dynamics a time step 
would correspond to a time increment At = 1/N [3D]. To 
speed up the arrival to the stationary state, we chose At 
larger, but still small enough so as to keep the process 
random. Typically, the probabilities involved in the evo- 
lution are such that max(£ At, pAt, gAt, hAt) < 0.1. 
We have checked that smaller probabilities do not affect 
the properties of the steady state, but only slow down the 
transient to reach it. We can recover the correct theoreti- 
cal expressions from our HMF analysis just by substitut- 
ing the model parameters by the rescaled values h/£, g/£, 
p/£, and, correspondingly, 77 = 1 + £ /p. 

6.1 Homogeneous networks: The Watts-Strogatz 
model 

As an example of a homogeneous network, we consider the 
small-world model proposed by Watts and Strogatz (WS) 
[23] . Networks in this model are generated as follows: The 
starting point is a ring with TV vertices, in which every 
vertex is symmetrically connected to its 2K nearest neigh- 
bors. Then, for every vertex, each edge connected to a 
clockwise neighbor is rewired to a randomly chosen vertex 
with probability p rw , and kept with probability 1 — p r w 
This procedure generates a graph with a degree distribu- 
tion that decays faster than exponentially for large k, and 
average degree (k) = 2K. We considered here WS net- 
works with p rw =l and K — 3. 

In Fig. [3] we show the numerical results obtained in 
WS homogeneous networks of size N — 10 5 , with fixed 



Fig. 3. Steady state fire density in the FFM on homogeneous 
WS networks of size N = 10 5 . (a) Raw data for different values 
of 77 as a function of g/£, and fixed h/i = 0.4. (b) Data collapse 
of the previous plots as given by Eq. (|47|l . The functions agree 
very well with the prediction by HFM theory for large values 
of the ratio g = r/g /£. The separation of the curves at low g 
indicates a possible non-mean-field regime. 

parameters £ = 10~ 3 and h/£ = 0.4, and varying values 
of p and g. We chose in particular a large value of h/l, 
larger than h c /£ = (k)^ 1 = 1/6, in order to avoid possi- 
ble problems in the vicinity of the critical point for small 
values of g. In this case, the theoretical prediction for the 
fire density is given by Eq. (]18|) . with the correct rescaling 
of parameters 

f _ 1 _ 9 + / 24 | + (1 + £ + Qh)2 

P F = - jt " — . (47) 

indicating, MS Ml' gued in Sec. [U that r\p F should be a scal- 
ing function of r\gj£. As we can observe in Fig. [3](b), the 
theoretical prediction is very well satisfied by numerical 
data for large values of rjg/i, collapsing all plots on the 
functional form of Eq. (|47|) when the appropriate rescal- 
ing is performed. However, for values of r\gjl < 10, we 
also observe a noticeable departure from the mean-field 
prediction, which is more conspicuous for large values of 
77. One could naively attribute this departure to a simple 
numerical artifact: since we have in general that p F ~ T)~ l , 
one could argue that, for fixed h, g, and £, larger values 
of rj lead to ever smaller fire densities. Therefore, for suf- 
ficiently large 77, we could expect such small fire density 
that its steady state determination will incur in numerical 
resolution problems, unless extremely large network sizes 
are considered. But here the minimal fire density wc find 
is of the order ~ 0.04 jj = 10~ 5 and so the number of 
burning trees still represents a significant fraction of the 
total population. 

A more thorough analysis, however, reveals that the 
actual reason of this departure is the failure of the mean- 
field assumption of lack of dynamical correlations between 
vertices discussed in Sec. [3] in the limit of small r/g/£ , and 
for large 77 (small p/ £)■ In this case, the fire density is very 
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Fig. 4. Dependence of a vertex's state on his neighbors' states 
in WS networks of size iV = 10 5 . The three graphs correspond 
to the three different curves presented in Fig(3] a) r\ = 1.1, 
b) i) = 2, c) X] = 11. The red curve is the fire density p F in 
each case. In the absence of state correlations, we should have 
P(F\E) = P{F\T) = P(F\F) = p F . 



small, and therefore burning vertices are very likely near- 
est neighbors of other fires, precisely those that originated 
them. This fact introduces correlations between the state 
of nearest neighbors which invalidate the whole mean-field 
approximation. We can check this argument in a homoge- 
neous network by comparing the probability P{a) that 
any vertex is in state a, with the conditional probability 
P(a\/3) that a vertex is in state a, provided it is nearest 
neighbor of a vertex in state (3. In Fig. [4] we compare the 
numerical values of P(F) = p F with the corresponding 
conditional probabilities P(F\E), P(F\T) and P(F\F) in 
simulations performed for different values of rj, plotted as 
a function of g/£. The figure shows that P(F\F) is clearly 
larger that p F in the areas of dissension: for small g/£, 
and especially for the largest values of 77. This results con- 
firms the presence of strong dynamical correlations be- 
tween vertices, and hints towards the failure of the HMF 
approximation in this region of the parameter space. 



6.2 Heterogeneous networks: The Barabasi-Albert 
model 

The Barabasi-Albert (BA) model is an algorithm to gen- 
erate growing SF networks with degree exponent 7=1, 
based on the preferential attachment paradigm [3]. This 
model is defined as follows: we start from a small num- 
ber m,Q of vertices, and at each time step, a new vertex is 
introduced, with m edges that are connected to old ver- 
tices i with probability PI (ki) = ki/^2jkj, where ki is 
the degree of the i th vertex. After iterating this procedure 
a large number of times, we obtain a network composed 
by N vertices, minimum degree m, fixed average degree 
(k) = 2m, degree distribution P(k) = 2m 2 k~ 3 and almost 
vanishing degree correlations [44 45j . The simulations con- 
sidered here were perform with m = 2 and a network size 
N = 10 6 . 



In the case of the BA network, we can check the ac- 
curacy of the exact HMF solution for 7 = 1, as given by 
Eqs. ([29]) and (f32|) . This is a nontrivial function of two 
variables, g and h. Therefore, to check it in a simple way, 
we focused on the small g and h regimes, in which expres- 
sions (|33|) and (|34[) should provide a good approximation. 
In particular, we computed the quantities dp F /dg and 
dp F /dh which, for small h/l and g can be approximated 

by 



df_ 

dg 



dh 



= 3 



8 - 10/im 



3=0 



h 2 m 2 ye.hm — 1 
8 - 2hm 



h 2 T 



(«*-!)■ 



hm ( e h m — 1 
2gm 



h=0 



(1+.9) 3 



(48) 



To obtain those quantities numerically, we used com- 
puter simulations to find the stationary fire density for 
several values of the variable with respect to which we 
took the partial derivative. Those values were chosen suf- 
ficiently small (g < 1CT 3 and h/l < 2 • 10~ 3 ) so that the 
expected second order term is negligible. We then checked 
that the points formed the expected straight line, on which 
we measured the slope. The error on this slope comes from 
statistical uncertainties on the measured fire density for 
the points considered. 

Figs. [5] and [6] show these quantities, computed from 
numerical simulations of the FFM in BA networks with 
I = 10~ 2 . In order to ensure that we are within the re- 
gion of validity of the HMF prediction, we chose a small 
value of 77 = 1.1 (jp/t = 10). The good agreement observed 
between numerical simulations and the theoretical HMF 
predictions in Eqs. (|48|) confirms the validity of the HMF 
analysis in this parameter regime. 



6.3 Heterogeneous networks: The uncorrelated 
configuration model 

To generate SF networks with an arbitrary degree ex- 
ponent 7 7^ 1, we used the uncorrelated configuration 
model (UCM) [46] . This model is defined as follows: We 
start from N initially disconnected vertices. Each vertex 
i is assigned a degree fcj, extracted from the probabil- 
ity distribution P(k) ~ fc~ 2 ~ 7 , subject to the constraints 
rn < ki < N 1 / 2 and £V k i even. Finally, the actual net- 
work is constructed by randomly connecting the vertices 
with Y]. kj/2 edges, respecting the preassigned degrees 
and avoiding multiple and self-connections. Using this al- 
gorithm, it is possible to create SF networks whose average 
maximum degree (or cut-off) scales as k c (N) ~ N 1 / 2 for 
any degree exponent 7, and which are completely uncorre- 
lated [47] . In the present simulations we chose a minimum 
degree m — 2, I = 10~ 2 and sizes up to = 10 6 . 
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a) 



b) 



Fig. 5. Slope ^j- as a function of h/l evaluated at g < 0.001 
on BA networks. The thick line is the mean field prediction for 
5 = 0, Eq.flM). 




Fig. 6. Slope -J- as a function of g evaluated at h/l < 0.002 
on BA networks. The thick line is the mean field prediction for 
h = 0, Eq.lpll. 
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Fig. 7. Steady state fire density in the FFM on heterogeneous 
UCM networks with 7=1/2 and size N = 10 6 . (a) Raw data 
for different values of 77 as a function oig/l, and fixed h/l = 0.4. 
(b) Data collapse of the previous plots as given in Eq.Jl4]). The 
full line is a numerical solution to Eq. ([26} plugged into (1291) 
with the appropriate parameters (m = 2, 7 = 1/2, h/l = 0.4). 
For small g = r]g/£, deviations from the mean field theory are 
observed. 
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In FigJT] we check the scaling of the fire density r\p F 
as a function of r/g/£ in UCM networks with degree ex- 
ponent 7 = 1/2. As in the case of the homogeneous WS 
networks, the data collapse is very good for large values of 
rig /I, fitting perfectly the theoretical prediction (full line) 
obtained by a numerical resolution of Eqs. |26|) and (|29|) . 
Again, deviations from the HMF prediction are observed 
for small r]g/£, which must be attributed to the presence 
of strong dynamical correlations between vertices, which 
invalidate the HMF approximation. 

The SF nature of this network model allows to explore 
the role of the degree in the establishment of dynamical 
correlations at small values of r\g/l. In order to do so, we 
concentrate in this case on the conditional probabilities 
Pk{a\f3) that a vertex of degree k be in state a, provided 
it is nearest neighbor of a vertex in state (3. In absence 
of dynamical correlations, we should expect Pk(a\/3) = p£ 
in the steady-state of the dynamics. In Fig. [8] we show 



Fig. 8. Dependence of a vertex's state on his neighbors' states, 
as a function of the degree in UCM networks with 7=5 an d 
size N = 10 6 . The thick line corresponds to the theoretical 
prediction , computed numerically from Eqs. (|26p and (|12p . 
Data for h/l = 0.4, g = 0.01. 

the theoretical value p*[ for g = 10~ 2 , calculated from 
Eqs. ((26|) and (fl"2|) . compared with the conditional rescaled 
probability T]Pk(F\F), evaluated from numerical simula- 
tions and plotted as a function of k for different values 
of 77. In the absence of dynamical correlations, we should 
observe the plots of ijPk{F\F) to collapse onto the theo- 
retical curve p F for the different values of rj. While this 
scenario is correct for 77 < 2, we observe very strong devi- 
ations from the HMF prediction for large 77, signaling the 
non-mean-ficld behavior of the FFM in SF networks. In 
particular, the conditional probability Pk(F\F) turns out 
to be larger than the HMF average value p F = rj^ 1 p F , 
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the difference increasing for small degree values. The dis- 
crepancy between Pk(F\F) and the average fire density, 
also observed in homogeneous networks (see Fig. [8]), is 
again due to the clustering at low fire densities of burning 
vertices in connected regions of the network. 

7 Conclusion 

In this paper we have presented a detailed analytical study 
of the forest fire model (FFM) in complex networks. From 
the perspective of the modelization of epidemic spreading 
processes, the FFM represents a generalization of several 
well-known epidemic models previously studied, which can 
be captured within the formalism of the FFM by the ap- 
propriate selection of representative parameters. Applying 
the now established HMF theory formalism, we have de- 
rived a set of rate equations in continuous time that rep- 
resent the dynamics of this model. Focusing in the long 
term steady state behavior, we have defined a set of alge- 
braic equations, whose analysis allows to discuss the role 
of the different parameters in the model. Thus, the rate at 
which empty sites become trees (the recovery rate of in- 
fected individuals) p, turns out to be absorbed in a rescal- 
ing of the fire density (density of infected individuals) and 
of the spontaneous ignition (spontaneous infection) rate, 
yielding in this way a fire density that is inversely propor- 
tional to l/p. In the case of homogeneous networks, we 
recover the results previously obtained in the mean-field 
analysis of the FFM [39] . In the case of heterogeneous SF 
networks, on the other hand, we have been able to provide 
exact explicit expressions for the fire density for a degree 
exponent 7 = 1, and approximate expressions for 7 < 1, 
valid for the cases of h or g very small. 

A comparison of these theoretical predictions with large 
scale simulations in homogeneous and heterogeneous net- 
works shows the HMF theory to provide a correct descrip- 
tion of the steady state of the FFM for large p (small 77) 
and g, a regime in which the average fire density is suffi- 
ciently large. For small p (large 77) and g, on the other 
hand, numerical simulations indicate the breakdown of 
HMF theory. The origin of this failure can be traced back 
to the build up of dynamical correlations between nearest 
neighbor vertices, correlations which in fact are expected 
to appear, due to the fact that fires accumulate with large 
probability in connected clusters, and are therefore not 
homogeneously distributed over the network as assumed 
by mean-field approaches. 

Apart from providing new insights into the behavior of 
a dynamical model relevant in epidemiological modeling, 
our results indicate a possible path to the understanding 
of the failure of HMF theory observed in other kinds of 
non-equilibrium processes in complex networks [TJ. 
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tional financial support through ICREA Academia, funded by 
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making available to us the Bionics cluster, and C. Castellano 
for helpful comments and discussions. 




Fig. 9. General dynamical process with a rest state A\. 



Appendix A: Rest states play no role in steady- 
state solutions 

Let us consider a general dynamical process taking place 
on 5 states where A\ is what we will call here a "rest 
state" : every site Ai can evolve to A\ only with some fixed 
rate £i and A\ evolves to A 2 with some other fixed rate p, 
see Fig. The dynamical equations of this process take 
the general form: 

'oi(i) = ^ i=2 hai{t) -pai(t) 
02 (t) = pai(t) + / a [oa(t), ...,os(*)] - ^o 2 (f) 

< Oj(t) =fi[a2(t),...,as(t)]-£iai(t) , (49) 

as(t) = f s [a 2 (t), ...,a s (t)) - £ s a s {t) 

> 1 = Ef=i o* 

for some functions /j [02(f), 05(f)], where Oj(f) is the 
probability of finding the system in state i at time f . Try- 
ing to solve this system for the stationary state by setting 
di(f) = Vi, one can get rid of the ai variable by sub- 
stituting its value from the first equation into the other 
ones, namely a\ = '^2 i £i/p ai- Then the last 5 — 1 equa- 
tions form a set of 5 — 1 independent equations with only 
5—1 unknowns: 

= Ef=2^ a ? + hi ®, -jOs) - ^a 2 (f) 

= fi(a 2 , a s ) - ti<H(t) ^ 

= /s(a2,...,as) -las 

Now if the introduced functions can be written as 
/ 4 (o 2 ,...,o s )= Yl <V " ' • "2 ■■■■■ "s ■ (51) 

ri2,—,ns 
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then, performing the change of variable 

<•, = VM, cT-' ns = <■; " II'//' • U = UvV 1 

(52) 

keeps all the first equalities in (|50|) true for any value of 
the constants In particular, choosing 



Tji = 1 + ii/p, 



(53) 



then all equalities arc verified and the last one becomes 
Ei=2 o* = 1- 

In other words, the new variables ai, together with 
Eqs. (|50[) above, describe the stationary solution of a dy- 
namical system with analogous dynamics as the first one 
([4"9")h up to a constant for some parameters, and without 
the rest state. Finding the stationary solution of the orig- 
inal S'-states system is thus equivalent to finding the one 
for the new system with 5—1 states. 

The new model has one parameter less, p, which has 
been absorbed in the redefinition of all other parameters 
and variables of the new system. So the stationary state 
of the first model is somehow independent on p. This can 
be stated as: 



'/-''- {{Vj , II C T " ' "" S ' }j ) 



Oi({tij""" nB ,lj}j) Indep. of pVi 



(54) 



Moreover, an upper bound for all densities a,; can be 
deduced directly from this argument: Since a, is the prob- 
ability to be in state i in the new system, then 



at < 1 =>■ a; < r)l 



(55) 



Notice that in the main text we chose the time scale 
so that i = 1. This has the effect of applying the ij factor 
to different coefficients (c.f. Eq. (fT0|) ). 



Appendix B: An explicit solution to a tran- 
scendental equation 



Let's consider the equation 

k — a + btt 1 



(56) 



for k, where a, b E M + are given and < 7 < 1 is 
fixed. This equation was already considered back in 1772 
by Lambert [38] . 

If a 7^ we can introduce the variable v = k/ 'a and 
see that this new variable depends on a unique variable 
p = bo*- 1 : 

i^l + Pv" 1 . (57) 

Similarly, if b 7^ 0, the new variable u = Kb"'- 1 only de- 

1 1 
pends on a = ab'- 1 = (it- 1 : 



u = a + u' 



(58) 



So if we know how to express explicitly v{(3) or u(a), 
we can conversely do so with K,(a,b). To achieve this, we 
use the Lagrange inversion theorem [49] which states that 
the reciprocal function x{y) of y — f(x) can be expressed 
around yo = f(xo) as 



x 



(59) 



with c„ = 



£,if /'(so) 9^0. 



dx" 1 \f(x)—y 0/ 

(1) Case a^O 

In this case we have (3{v) = v^ 7 (v — 1) and so around 
vq = 1 and /3q = 0, /3'(1) 7^ and we can calculate: 



v - v 



dv 11 - 1 \(3(v) ~ fa 

r(n 7 + 1) 
r(n(7-l) + 2)' 



Thus we can write 



r(n7 + i) 



^ i r(n( 7 -l) + 2) J T(n+l) 



(60) 
(61) 

(62) 



This formulation of the function v(/3) is true as long as 
the series converges. To know what interval it corresponds 
to, we calculate the convergence radius f3 c following its 
definition: 




sin(ri7r(l — 7)) 



(63) 

r(n 7 + !)/>(!- 7) -1) 



r(n + l) 



= lim 



B /r(n7+l)r(n(l-7)-l) 



r(n + l) 



(64) 
(65) 



where we have used the Euler reflection formula to avoid 
an undefined valued of the Gamma function [50]. Using 
the Stirling formula T(z) = ^r(f) Z (l + 0{\)) and 
liirin—xjo \fn = 1 we find 



Pc 1 =7 7 (l-7) 



1-7 



(66) 



So Eq. ([51]) can be used to express k explicitly, but not 
for any value of its arguments. 
(2) Case 6^0: 

Let's find now a series representation of function u(a). 
To simplify the calculation we introduce in this case w = 
it 1-7 so that 

1 7 
w x --t = a + w 1 —' . (67) 

We can now directly apply the same steps as before to 
find a power series expression for the reciprocal function 
of a(w) = w 1 ~'< (w — 1) around wo = 1, ot = 0, since 
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a'(l)^0: 



W — Wq 



dw n 1 \a(w) — Qo 



V 1— 7 



(68) 
(69) 



Thus we have 



u{a) = 



7 



i) 



^r(-i|2)r(,ii) 



(70) 



9. 
10. 

11. 

12. 

13. 

14. 

15. 
16. 



Again we can calculate the convergence radius and see ^' • 
_j_ i 18. 

that it is a c = (1 — 7)7 1 ~~> = /3 C T . In other words, 

we found a complete description of the solution re(<x, 6) of 19 

Eq. (|56|) : depending on whether the value of (3 is smaller or 

greater than /3 C , development (|62|) or ([70|l converges and 20. 

gives an explicit value for n(a,b): 



E 

n=0 



r(»7+i) 



r(n( 7 -l)+2)r(n+l) 



if /? < (3 C 



re(a, 6) = < 



with 



if /3 > /3 C 

/? c = 7^(i-7)" (1 " 7) e]l,2[. 



(71) 
(72) 



21. 

22. 

23. 
24. 

25. 
26. 
27. 
28. 
29. 
30. 



Since k depends smoothly on a and b, as can be seen from 31 
Eq. (|56p . the case (3 = (3 C is deducible as the limit of any 
of the two series. 32 
This expression gives the value of the function n(a,b) ^3 
whenever either a or b is different from zero. The solution 
in the case a = b = is trivially x = 0. 



References 



1. A. Barrat, M. Barthelemy, A. Vespignani, Dynamical Pro 
cesses on Complex Networks (Cambridge University Press, 
Cambridge, 2008) 

2. S.N. Dorogovtsev, A.V. Goltsev, J.F.F. Mendes, Rev. 
Mod. Phys. 80, 1275 (2008) 37. 

3. A.L. Barabasi, R. Albert, Science 286, 509 (1999) 

4. G. Caldarelli, Scale-Free Networks: Complex Webs in Na- 38. 
ture and Technology (Oxford University Press, Oxford, 
2007) 39. 

5. D.S. Callaway, M.E.J. Newman, S.H. Strogatz, D.J. Watts, 
Phys. Rev. Lett. 85, 5468 (2000) 40. 

6. R. Cohen, K. Erez, D. ben Avraham, S. Havlin, Phys. Rev. 
Lett. 86, 3682 (2001) 

7. C. Castellano, R. Pastor-Satorras, Phys. Rev. Lett. 96, 41. 
038701 (2006) 

8. K.I. Goh, D.S. Lee, B. Kahng, D. Kim, Phys. Rev. Lett. 42. 
91, 148701 (2003) 



Y. Moreno, A. Vazquez, Europhys. Lett. 57, 765 (2002) 
A. Arenas, A. Dfaz-Guilera, J. Kurths, Y. Moreno, 
C. Zhou, Phys. Rep. 469, 93 (2008) 

C. Castellano, S. Fortunato, V. Loreto, Rev. Mod. Phys. 
81, 591 (2008) 

M.J. Keeling, K.T.D. Eames, J. R. Soc. Interface 2, 295 
(2005) 

F. Liljeros, C.R. Edling, L.A.N. Amaral, H.E. Stanley, 
Y. Aberg, Nature 411, 907 (2001) 

R. Pastor-Satorras, A. Vespignani, Phys. Rev. Lett. 
86(14), 3200 (2001) 

A. L. Lloyd, R.M. May, Science 292, 1316 (2001) 

Y. Moreno, R. Pastor-Satorras, A. Vespignani, Eur. Phys. 
J. B 26, 521 (2002) 

M.E.J. Newman, Phys. Rev. E 66, 016128 (2002) 

R. Pastor-Satorras, A. Vespignani, Phys. Rev. E 65, 

036104 (2001) 

R. Cohen, S. Havlin, D. ben- Avraham, Physical Review 
Letters 91, 247901 (2003) 

J.G. Gardenes, V. Latora, Y. Moreno, E. Profumo, Proc. 

Natl. Acad. Sci. USA 105, 1399 (2008) 

V. Colizza, A. Barrat, M. Barthelemy, A. Vespignani, Proc. 

Natl. Acad. Sci. USA 103, 2015 (2006) 

R.M. Anderson, R.M. May, Infectious diseases in humans 

(Oxford University Press, Oxford, 1992) 

D. J. Watts, S.H. Strogatz, Nature 393, 440 (1998) 

M. Boguna, C. Castellano, R. Pastor-Satorras, Phys. Rev. 
E 79, 036110 (2009) 

P. Bak, K. Chen, C. Tang, Phys. Lett. A 147, 297 (1992) 

B. Drossel, F. Schwabl, Phys. Rev. Lett. 69, 1629 (1992) 
P. Grassberger, New Journal of Physics 4, 17 (2002) 
J.A. Bonachela, M.A. Munoz, J STAT p. P09009 (2009) 

C. Rhodes, R. Anderson, Nature 381, 600 (1996) 

G. Abramson, M. Kuperman, Phys. Rev. Lett. 86, 2909 
(2001) 

CM. M. Miiller-Linow, M.T. Hiitt, Phys. Rev. E 74, 
016112 (2006) 

J. Liu, Y. Tang, Z. Yang, J. Stat. Mech. p. P08008 (2004) 
S. Peng, Y. Li, B. Zheng, Steady states and critical behav- 
ior of epidemic spreading on complex networks (WCICA 
2008. 7th World Congress on Intelligent Control and Au- 
tomation, 2008), pp. 3481-3486 

R. Pastor-Satorras, A. Vazquez, A. Vespignani, Phys. Rev. 
Lett. 87, 258701 (2001) 

H. J. Jensen, Self- Organized Criticality (Cambridge Uni- 
versity Press, Cambridge, 1998) 

M. Boguna, R. Pastor-Satorras, A. Vespignani, in Statisti- 
cal Mechanics of Complex Networks, edited by R. Pastor- 
Satorras, J.M. Rubf, A. Dfaz-Guilera (Springer Verlag, 
Berlin, 2003), Vol. 625 of Lecture Notes in Physics 
R. Pastor-Satorras, A. Vespignani, Phys. Rev. E 63, 
066117 (2001) 

M. Catanzaro, M. Boguna, R. Pastor-Satorras, Phys. Rev. 
E 71, 056104 (2005) 

K. Christensen, H. Flyvbjerg, Z. Olami, Phys. Rev. Lett. 
71, 2737 (1993) 

J. Marro, R. Dickman, Nonequilibrium phase transitions 
in lattice models (Cambridge University Press, Cambridge, 
1999) 

M. Boguna, R. Pastor-Satorras, Phys. Rev. E 66, 047104 
(2002) 

M. Abramowitz, LA. Stegun, Handbook of mathematical 
functions. (Dover, New York, 1972) 



J.-D. Bancal and R. Pastor-Satorras: Steady-State Dynamics of the FFM on Complex Networks 



43. R.M. Corless, G.H. Gonnet, D.E.G. Hare, D.J. Jeffrey, 
D.E. Knuth, Advances in Computational Mathematics 5, 
329 (1996) 

44. A. Vazquez, R. Pastor-Satorras, A. Vespignani, Phys. Rev. 
E 65, 066130 (2002) 

45. A. Barrat, R. Pastor-Satorras, Phys. Rev. E 71, 036127 
(2005) 

46. M. Catanzaro, M. Boguna, R. Pastor-Satorras, Phys. Rev. 
E 71, 027103 (2005) 

47. M. Boguna, R. Pastor-Satorras, A. Vespignani, Euro. 
Phys. J. B 38, 205 (2004) 

48. J.H. Lambert, Nouveaux memoires de l'Academie royale 
des sciences et belles-lettres 1 (1772) 

49. E. Goursat, Course in Mathematical Analysis, Vol. 2: 
Functions of a Complex Variable & Differential Equations 
(Dover, New York, 1959) 

50. J. Havil, Exploring Euler's Constant (Princeton University 
Press, Princeton, NJ, 2003) 



